// Copyright (C) 2000, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#if defined(_MSC_VER)
// Turn off compiler warning about long names
#pragma warning(disable : 4786)
#endif

#include <cassert>
#include <cstdio>

#include "CoinTypes.h"
#include "CoinFloatEqual.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinIndexedVector.hpp"

//#############################################################################
#define WARN_USELESS 0
void CoinIndexedVector::clear()
{
#if WARN_USELESS
  int nNonZero = 0;
  for (int i = 0; i < capacity_; i++) {
    if (elements_[i])
      nNonZero++;
  }
  assert(nNonZero <= nElements_);
#if WARN_USELESS > 1
  if (nNonZero != nElements_)
    printf("Vector said it had %d nonzeros - only had %d\n",
      nElements_, nNonZero);
#endif
  if (!nNonZero && nElements_)
    printf("Vector said it had %d nonzeros - but is already empty\n",
      nElements_);
#endif
  assert(nElements_ <= capacity_);
  if (!packedMode_) {
#ifndef NDEBUG
    for (int i = 0; i < nElements_; i++)
      assert(indices_[i] >= 0 && indices_[i] < capacity_);
#endif
    if (3 * nElements_ < capacity_) {
      int i = 0;
      if ((nElements_ & 1) != 0) {
        elements_[indices_[0]] = 0.0;
        i = 1;
      }
      for (; i < nElements_; i += 2) {
        int i0 = indices_[i];
        int i1 = indices_[i + 1];
        elements_[i0] = 0.0;
        elements_[i1] = 0.0;
      }
    } else {
      CoinZeroN(elements_, capacity_);
    }
  } else {
    CoinZeroN(elements_, nElements_);
  }
  nElements_ = 0;
  packedMode_ = false;
  //checkClear();
}
void CoinIndexedVector::reallyClear()
{
  CoinZeroN(elements_, capacity_);
  nElements_ = 0;
  packedMode_ = false;
}
//#############################################################################

void CoinIndexedVector::empty()
{
  delete[] indices_;
  indices_ = NULL;
  if (elements_)
    delete[](elements_ - offset_);
  elements_ = NULL;
  nElements_ = 0;
  capacity_ = 0;
  packedMode_ = false;
}

//#############################################################################

CoinIndexedVector &
CoinIndexedVector::operator=(const CoinIndexedVector &rhs)
{
  if (this != &rhs) {
    clear();
    packedMode_ = rhs.packedMode_;
    if (!packedMode_)
      gutsOfSetVector(rhs.capacity_, rhs.nElements_,
        rhs.indices_, rhs.elements_);
    else
      gutsOfSetPackedVector(rhs.capacity_, rhs.nElements_,
        rhs.indices_, rhs.elements_);
  }
  return *this;
}
/* Copy the contents of one vector into another.  If multiplier is 1
   It is the equivalent of = but if vectors are same size does
   not re-allocate memory just clears and copies */
void CoinIndexedVector::copy(const CoinIndexedVector &rhs, double multiplier)
{
  if (capacity_ == rhs.capacity_) {
    // can do fast
    clear();
    packedMode_ = rhs.packedMode_;
    nElements_ = 0;
    if (!packedMode_) {
      for (int i = 0; i < rhs.nElements_; i++) {
        int index = rhs.indices_[i];
        double value = rhs.elements_[index] * multiplier;
        if (fabs(value) < COIN_INDEXED_TINY_ELEMENT)
          value = COIN_INDEXED_REALLY_TINY_ELEMENT;
        elements_[index] = value;
        indices_[nElements_++] = index;
      }
    } else {
      for (int i = 0; i < rhs.nElements_; i++) {
        int index = rhs.indices_[i];
        double value = rhs.elements_[i] * multiplier;
        if (fabs(value) < COIN_INDEXED_TINY_ELEMENT)
          value = COIN_INDEXED_REALLY_TINY_ELEMENT;
        elements_[nElements_] = value;
        indices_[nElements_++] = index;
      }
    }
  } else {
    // do as two operations
    *this = rhs;
    (*this) *= multiplier;
  }
}

//#############################################################################
#ifndef CLP_NO_VECTOR
CoinIndexedVector &
CoinIndexedVector::operator=(const CoinPackedVectorBase &rhs)
{
  clear();
  packedMode_ = false;
  gutsOfSetVector(rhs.getNumElements(),
    rhs.getIndices(), rhs.getElements());
  return *this;
}
#endif
//#############################################################################

void CoinIndexedVector::borrowVector(int size, int numberIndices, int *inds, double *elems)
{
  empty();
  capacity_ = size;
  nElements_ = numberIndices;
  indices_ = inds;
  elements_ = elems;

  // whole point about borrowvector is that it is lightweight so no testing is done
}

//#############################################################################

void CoinIndexedVector::returnVector()
{
  indices_ = NULL;
  elements_ = NULL;
  nElements_ = 0;
  capacity_ = 0;
  packedMode_ = false;
}

//#############################################################################

void CoinIndexedVector::setVector(int size, const int *inds, const double *elems)
{
  clear();
  gutsOfSetVector(size, inds, elems);
}
//#############################################################################

void CoinIndexedVector::setVector(int size, int numberIndices, const int *inds, const double *elems)
{
  clear();
  gutsOfSetVector(size, numberIndices, inds, elems);
}
//#############################################################################

void CoinIndexedVector::setConstant(int size, const int *inds, double value)
{
  clear();
  gutsOfSetConstant(size, inds, value);
}

//#############################################################################

void CoinIndexedVector::setFull(int size, const double *elems)
{
  // Clear out any values presently stored
  clear();

#ifndef COIN_FAST_CODE
  if (size < 0)
    throw CoinError("negative number of indices", "setFull", "CoinIndexedVector");
#endif
  reserve(size);
  nElements_ = 0;
  // elements_ array is all zero
  int i;
  for (i = 0; i < size; i++) {
    int indexValue = i;
    if (fabs(elems[i]) >= COIN_INDEXED_TINY_ELEMENT) {
      elements_[indexValue] = elems[i];
      indices_[nElements_++] = indexValue;
    }
  }
}
//#############################################################################

/** Access the i'th element of the full storage vector.  */
double &
  CoinIndexedVector::operator[](int index) const
{
  assert(!packedMode_);
#ifndef COIN_FAST_CODE
  if (index >= capacity_)
    throw CoinError("index >= capacity()", "[]", "CoinIndexedVector");
  if (index < 0)
    throw CoinError("index < 0", "[]", "CoinIndexedVector");
#endif
  double *where = elements_ + index;
  return *where;
}
//#############################################################################

void CoinIndexedVector::setElement(int index, double element)
{
#ifndef COIN_FAST_CODE
  if (index >= nElements_)
    throw CoinError("index >= size()", "setElement", "CoinIndexedVector");
  if (index < 0)
    throw CoinError("index < 0", "setElement", "CoinIndexedVector");
#endif
  elements_[indices_[index]] = element;
}

//#############################################################################

void CoinIndexedVector::insert(int index, double element)
{
#ifndef COIN_FAST_CODE
  if (index < 0)
    throw CoinError("index < 0", "setElement", "CoinIndexedVector");
#endif
  if (index >= capacity_)
    reserve(index + 1);
#ifndef COIN_FAST_CODE
  if (elements_[index])
    throw CoinError("Index already exists", "insert", "CoinIndexedVector");
#endif
  indices_[nElements_++] = index;
  elements_[index] = element;
}

//#############################################################################

void CoinIndexedVector::add(int index, double element)
{
#ifndef COIN_FAST_CODE
  if (index < 0)
    throw CoinError("index < 0", "setElement", "CoinIndexedVector");
#endif
  if (index >= capacity_)
    reserve(index + 1);
  if (elements_[index]) {
    element += elements_[index];
    if (fabs(element) >= COIN_INDEXED_TINY_ELEMENT) {
      elements_[index] = element;
    } else {
      elements_[index] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    }
  } else if (fabs(element) >= COIN_INDEXED_TINY_ELEMENT) {
    indices_[nElements_++] = index;
    assert(nElements_ <= capacity_);
    elements_[index] = element;
  }
}

//#############################################################################

int CoinIndexedVector::clean(double tolerance)
{
  int number = nElements_;
  int i;
  nElements_ = 0;
  assert(!packedMode_);
  for (i = 0; i < number; i++) {
    int indexValue = indices_[i];
    if (fabs(elements_[indexValue]) >= tolerance) {
      indices_[nElements_++] = indexValue;
    } else {
      elements_[indexValue] = 0.0;
    }
  }
  return nElements_;
}
#ifndef NDEBUG
//#############################################################################
// For debug check vector is clear i.e. no elements
void CoinIndexedVector::checkClear()
{
#ifndef NDEBUG
  //printf("checkClear %p\n",this);
  assert(!nElements_);
  //assert(!packedMode_);
  int i;
  for (i = 0; i < capacity_; i++) {
    assert(!elements_[i]);
  }
  // check mark array zeroed
  char *mark = reinterpret_cast< char * >(indices_ + capacity_);
  for (i = 0; i < capacity_; i++) {
    assert(!mark[i]);
  }
#else
  if (nElements_) {
    printf("%d nElements_ - checkClear\n", nElements_);
    abort();
  }
  if (packedMode_) {
    printf("packed mode when empty - checkClear\n");
    abort();
  }
  int i;
  int n = 0;
  int k = -1;
  for (i = 0; i < capacity_; i++) {
    if (elements_[i]) {
      n++;
      if (k < 0)
        k = i;
    }
  }
  if (n) {
    printf("%d elements, first %d - checkClear\n", n, k);
    abort();
  }
#endif
}
// For debug check vector is clean i.e. elements match indices
void CoinIndexedVector::checkClean()
{
  //printf("checkClean %p\n",this);
  int i;
  if (packedMode_) {
    for (i = 0; i < nElements_; i++)
      assert(elements_[i]);
    for (; i < capacity_; i++)
      assert(!elements_[i]);
  } else {
    double *copy = new double[capacity_];
    CoinMemcpyN(elements_, capacity_, copy);
    for (i = 0; i < nElements_; i++) {
      int indexValue = indices_[i];
      assert(copy[indexValue]);
      copy[indexValue] = 0.0;
    }
    for (i = 0; i < capacity_; i++)
      assert(!copy[i]);
    delete[] copy;
  }
#ifndef NDEBUG
  // check mark array zeroed
  char *mark = reinterpret_cast< char * >(indices_ + capacity_);
  for (i = 0; i < capacity_; i++) {
    assert(!mark[i]);
  }
#endif
}
#endif
//#############################################################################
#ifndef CLP_NO_VECTOR
void CoinIndexedVector::append(const CoinPackedVectorBase &caboose)
{
  const int cs = caboose.getNumElements();

  const int *cind = caboose.getIndices();
  const double *celem = caboose.getElements();
  int maxIndex = -1;
  int i;
  for (i = 0; i < cs; i++) {
    int indexValue = cind[i];
    if (indexValue < 0)
      throw CoinError("negative index", "append", "CoinIndexedVector");
    if (maxIndex < indexValue)
      maxIndex = indexValue;
  }
  reserve(maxIndex + 1);
  bool needClean = false;
  int numberDuplicates = 0;
  for (i = 0; i < cs; i++) {
    int indexValue = cind[i];
    if (elements_[indexValue]) {
      numberDuplicates++;
      elements_[indexValue] += celem[i];
      if (fabs(elements_[indexValue]) < COIN_INDEXED_TINY_ELEMENT)
        needClean = true; // need to go through again
    } else {
      if (fabs(celem[i]) >= COIN_INDEXED_TINY_ELEMENT) {
        elements_[indexValue] = celem[i];
        indices_[nElements_++] = indexValue;
      }
    }
  }
  if (needClean) {
    // go through again
    int size = nElements_;
    nElements_ = 0;
    for (i = 0; i < size; i++) {
      int indexValue = indices_[i];
      double value = elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        indices_[nElements_++] = indexValue;
      } else {
        elements_[indexValue] = 0.0;
      }
    }
  }
  if (numberDuplicates)
    throw CoinError("duplicate index", "append", "CoinIndexedVector");
}
#endif
//#############################################################################

void CoinIndexedVector::swap(int i, int j)
{
  if (i >= nElements_)
    throw CoinError("index i >= size()", "swap", "CoinIndexedVector");
  if (i < 0)
    throw CoinError("index i < 0", "swap", "CoinIndexedVector");
  if (j >= nElements_)
    throw CoinError("index j >= size()", "swap", "CoinIndexedVector");
  if (j < 0)
    throw CoinError("index j < 0", "swap", "CoinIndexedVector");

  // Swap positions i and j of the
  // indices array

  int isave = indices_[i];
  indices_[i] = indices_[j];
  indices_[j] = isave;
}

//#############################################################################

void CoinIndexedVector::truncate(int n)
{
  reserve(n);
}

//#############################################################################

void CoinIndexedVector::operator+=(double value)
{
  assert(!packedMode_);
  int i, indexValue;
  for (i = 0; i < nElements_; i++) {
    indexValue = indices_[i];
    double newValue = elements_[indexValue] + value;
    if (fabs(newValue) >= COIN_INDEXED_TINY_ELEMENT)
      elements_[indexValue] = newValue;
    else
      elements_[indexValue] = COIN_INDEXED_REALLY_TINY_ELEMENT;
  }
}

//-----------------------------------------------------------------------------

void CoinIndexedVector::operator-=(double value)
{
  assert(!packedMode_);
  int i, indexValue;
  for (i = 0; i < nElements_; i++) {
    indexValue = indices_[i];
    double newValue = elements_[indexValue] - value;
    if (fabs(newValue) >= COIN_INDEXED_TINY_ELEMENT)
      elements_[indexValue] = newValue;
    else
      elements_[indexValue] = COIN_INDEXED_REALLY_TINY_ELEMENT;
  }
}

//-----------------------------------------------------------------------------

void CoinIndexedVector::operator*=(double value)
{
  assert(!packedMode_);
  int i, indexValue;
  for (i = 0; i < nElements_; i++) {
    indexValue = indices_[i];
    double newValue = elements_[indexValue] * value;
    if (fabs(newValue) >= COIN_INDEXED_TINY_ELEMENT)
      elements_[indexValue] = newValue;
    else
      elements_[indexValue] = COIN_INDEXED_REALLY_TINY_ELEMENT;
  }
}

//-----------------------------------------------------------------------------

void CoinIndexedVector::operator/=(double value)
{
  assert(!packedMode_);
  int i, indexValue;
  for (i = 0; i < nElements_; i++) {
    indexValue = indices_[i];
    double newValue = elements_[indexValue] / value;
    if (fabs(newValue) >= COIN_INDEXED_TINY_ELEMENT)
      elements_[indexValue] = newValue;
    else
      elements_[indexValue] = COIN_INDEXED_REALLY_TINY_ELEMENT;
  }
}
//#############################################################################

void CoinIndexedVector::reserve(int n)
{
  assert (!packedMode_||!nElements_);
#if OLD_COININDEXED_RESERVE_CODE
  /* Code was modified in 2017 to fix a bug
     I can not see justification for the changes then, so going back 
     to approximately old way.  
     Leaving this in for another few years - just in case. */
  int i;
  int nPlus;
  if (sizeof(int) == 4 * sizeof(char))
    nPlus = (n + 3) >> 2;
  else
    nPlus = (n + 7) >> 4;
#ifdef COIN_AVX2
  nPlus += 16;
#endif
  // don't make allocated space smaller but do take off values
  if (n + nPlus < capacity_) {
#ifndef COIN_FAST_CODE
    if (n < 0)
      throw CoinError("negative capacity", "reserve", "CoinIndexedVector");
#endif
    int nNew = 0;
    for (i = 0; i < nElements_; i++) {
      int indexValue = indices_[i];
      if (indexValue < n) {
        indices_[nNew++] = indexValue;
      } else {
        elements_[indexValue] = 0.0;
      }
    }
    nElements_ = nNew;
  } else if (n > capacity_) {

    // save pointers to existing data
    int *tempIndices = indices_;
    double *tempElements = elements_;
    double *delTemp = elements_ - offset_;

    // allocate new space
    indices_ = new int[n + nPlus];
    CoinZeroN(indices_ + n, nPlus);
    // align on 64 byte boundary
    double *temp = new double[n + 9 + nPlus];
    offset_ = 0;
    CoinInt64 xx = reinterpret_cast< CoinInt64 >(temp);
    int iBottom = static_cast< int >(xx & 63);
    //if (iBottom)
    offset_ = (64 - iBottom) >> 3;
    elements_ = temp + offset_;
    ;

    // copy data to new space
    // and zero out part of array
    if (nElements_ > 0) {
      CoinMemcpyN(tempIndices, nElements_, indices_);
      CoinMemcpyN(tempElements, capacity_, elements_);
      CoinZeroN(elements_ + capacity_, n - capacity_);
    } else {
      CoinZeroN(elements_, n);
    }
    capacity_ = n;

    // free old data
    if (tempElements)
      delete[] delTemp;
    delete[] tempIndices;
  }
#else
  // extra space for char array off end and for alignment
  int nPlus;
  if (sizeof(int) == 4 * sizeof(char))
    nPlus = (n + 3) >> 2;
  else
    nPlus = (n + 7) >> 4;
#ifdef COIN_AVX2
  nPlus += 16;
#endif
  // don't make allocated space smaller but do take off values
  if (n < capacity_) {
#ifndef COIN_FAST_CODE
    if (n < 0)
      throw CoinError("negative capacity", "reserve", "CoinIndexedVector");
#endif
    int nNew = 0;
    for (int i = 0; i < nElements_; i++) {
      int indexValue = indices_[i];
      if (indexValue < n) {
        indices_[nNew++] = indexValue;
      } else {
        elements_[indexValue] = 0.0;
      }
    }
    nElements_ = nNew;
  } else if (n > capacity_) {
    // save pointers to existing data
    int *tempIndices = indices_;
    double *tempElements = elements_;
    // true address for delete []
    double *delTemp = elements_ - offset_;

    // allocate new space
    indices_ = new int[n + nPlus];
    // align doubles on 64 byte boundary
#ifndef COIN_AVX2
    double *temp = new double[n + 9];
#else
    double *temp = new double[n + 11]; // allow extra
#endif
    offset_ = 0;
    CoinInt64 xx = reinterpret_cast< CoinInt64 >(temp);
    int iBottom = static_cast< int >(xx & 63);
    offset_ = (64 - iBottom) >> 3;
    elements_ = temp + offset_;

    // copy data to new space
    // and zero out part of array
    if (nElements_ > 0) {
      CoinMemcpyN(tempIndices, nElements_, indices_);
      CoinMemcpyN(tempElements, capacity_, elements_);
      CoinZeroN(elements_ + capacity_, n - capacity_);
    } else {
      CoinZeroN(elements_, n);
    }
    // free old data
    delete[] delTemp;
    delete[] tempIndices;
    capacity_ = n;
  }
  // be on safe side - zero out hidden char array
  if (indices_)
    CoinZeroN(indices_ + n, nPlus);
#endif
}

//#############################################################################

CoinIndexedVector::CoinIndexedVector()
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
}

CoinIndexedVector::CoinIndexedVector(int size)
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
  // Get space
  reserve(size);
}

//-----------------------------------------------------------------------------

CoinIndexedVector::CoinIndexedVector(int size,
  const int *inds, const double *elems)
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
  gutsOfSetVector(size, inds, elems);
}

//-----------------------------------------------------------------------------

CoinIndexedVector::CoinIndexedVector(int size,
  const int *inds, double value)
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
  gutsOfSetConstant(size, inds, value);
}

//-----------------------------------------------------------------------------

CoinIndexedVector::CoinIndexedVector(int size, const double *element)
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
  setFull(size, element);
}

//-----------------------------------------------------------------------------
#ifndef CLP_NO_VECTOR
CoinIndexedVector::CoinIndexedVector(const CoinPackedVectorBase &rhs)
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
  gutsOfSetVector(rhs.getNumElements(),
    rhs.getIndices(), rhs.getElements());
}
#endif
//-----------------------------------------------------------------------------

CoinIndexedVector::CoinIndexedVector(const CoinIndexedVector &rhs)
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
  if (!rhs.packedMode_)
    gutsOfSetVector(rhs.capacity_, rhs.nElements_, rhs.indices_, rhs.elements_);
  else
    gutsOfSetPackedVector(rhs.capacity_, rhs.nElements_, rhs.indices_, rhs.elements_);
}

//-----------------------------------------------------------------------------

CoinIndexedVector::CoinIndexedVector(const CoinIndexedVector *rhs)
  : indices_(NULL)
  , elements_(NULL)
  , nElements_(0)
  , capacity_(0)
  , offset_(0)
  , packedMode_(false)
{
  if (!rhs->packedMode_)
    gutsOfSetVector(rhs->capacity_, rhs->nElements_, rhs->indices_, rhs->elements_);
  else
    gutsOfSetPackedVector(rhs->capacity_, rhs->nElements_, rhs->indices_, rhs->elements_);
}

//-----------------------------------------------------------------------------

CoinIndexedVector::~CoinIndexedVector()
{
  delete[] indices_;
  if (elements_)
    delete[](elements_ - offset_);
}
//#############################################################################
//#############################################################################

/// Return the sum of two indexed vectors
CoinIndexedVector
CoinIndexedVector::operator+(
  const CoinIndexedVector &op2)
{
  assert(!packedMode_);
  int i;
  int nElements = nElements_;
  int capacity = std::max(capacity_, op2.capacity_);
  CoinIndexedVector newOne(*this);
  newOne.reserve(capacity);
  bool needClean = false;
  // new one now can hold everything so just modify old and add new
  for (i = 0; i < op2.nElements_; i++) {
    int indexValue = op2.indices_[i];
    double value = op2.elements_[indexValue];
    double oldValue = elements_[indexValue];
    if (!oldValue) {
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        newOne.elements_[indexValue] = value;
        newOne.indices_[nElements++] = indexValue;
      }
    } else {
      value += oldValue;
      newOne.elements_[indexValue] = value;
      if (fabs(value) < COIN_INDEXED_TINY_ELEMENT) {
        needClean = true;
      }
    }
  }
  newOne.nElements_ = nElements;
  if (needClean) {
    // go through again
    newOne.nElements_ = 0;
    for (i = 0; i < nElements; i++) {
      int indexValue = newOne.indices_[i];
      double value = newOne.elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        newOne.indices_[newOne.nElements_++] = indexValue;
      } else {
        newOne.elements_[indexValue] = 0.0;
      }
    }
  }
  return newOne;
}

/// Return the difference of two indexed vectors
CoinIndexedVector
CoinIndexedVector::operator-(
  const CoinIndexedVector &op2)
{
  assert(!packedMode_);
  int i;
  int nElements = nElements_;
  int capacity = std::max(capacity_, op2.capacity_);
  CoinIndexedVector newOne(*this);
  newOne.reserve(capacity);
  bool needClean = false;
  // new one now can hold everything so just modify old and add new
  for (i = 0; i < op2.nElements_; i++) {
    int indexValue = op2.indices_[i];
    double value = op2.elements_[indexValue];
    double oldValue = elements_[indexValue];
    if (!oldValue) {
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        newOne.elements_[indexValue] = -value;
        newOne.indices_[nElements++] = indexValue;
      }
    } else {
      value = oldValue - value;
      newOne.elements_[indexValue] = value;
      if (fabs(value) < COIN_INDEXED_TINY_ELEMENT) {
        needClean = true;
      }
    }
  }
  newOne.nElements_ = nElements;
  if (needClean) {
    // go through again
    newOne.nElements_ = 0;
    for (i = 0; i < nElements; i++) {
      int indexValue = newOne.indices_[i];
      double value = newOne.elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        newOne.indices_[newOne.nElements_++] = indexValue;
      } else {
        newOne.elements_[indexValue] = 0.0;
      }
    }
  }
  return newOne;
}

/// Return the element-wise product of two indexed vectors
CoinIndexedVector
  CoinIndexedVector::operator*(
    const CoinIndexedVector &op2)
{
  assert(!packedMode_);
  int i;
  int nElements = nElements_;
  int capacity = std::max(capacity_, op2.capacity_);
  CoinIndexedVector newOne(*this);
  newOne.reserve(capacity);
  bool needClean = false;
  // new one now can hold everything so just modify old and add new
  for (i = 0; i < op2.nElements_; i++) {
    int indexValue = op2.indices_[i];
    double value = op2.elements_[indexValue];
    double oldValue = elements_[indexValue];
    if (oldValue) {
      value *= oldValue;
      newOne.elements_[indexValue] = value;
      if (fabs(value) < COIN_INDEXED_TINY_ELEMENT) {
        needClean = true;
      }
    }
  }

  newOne.nElements_ = nElements;

  if (needClean) {
    // go through again
    newOne.nElements_ = 0;
    for (i = 0; i < nElements; i++) {
      int indexValue = newOne.indices_[i];
      double value = newOne.elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        newOne.indices_[newOne.nElements_++] = indexValue;
      } else {
        newOne.elements_[indexValue] = 0.0;
      }
    }
  }
  return newOne;
}

/// Return the element-wise ratio of two indexed vectors
CoinIndexedVector
CoinIndexedVector::operator/(const CoinIndexedVector &op2)
{
  assert(!packedMode_);
  // I am treating 0.0/0.0 as 0.0
  int i;
  int nElements = nElements_;
  int capacity = std::max(capacity_, op2.capacity_);
  CoinIndexedVector newOne(*this);
  newOne.reserve(capacity);
  bool needClean = false;
  // new one now can hold everything so just modify old and add new
  for (i = 0; i < op2.nElements_; i++) {
    int indexValue = op2.indices_[i];
    double value = op2.elements_[indexValue];
    double oldValue = elements_[indexValue];
    if (oldValue) {
      if (!value)
        throw CoinError("zero divisor", "/", "CoinIndexedVector");
      value = oldValue / value;
      newOne.elements_[indexValue] = value;
      if (fabs(value) < COIN_INDEXED_TINY_ELEMENT) {
        needClean = true;
      }
    }
  }

  newOne.nElements_ = nElements;

  if (needClean) {
    // go through again
    newOne.nElements_ = 0;
    for (i = 0; i < nElements; i++) {
      int indexValue = newOne.indices_[i];
      double value = newOne.elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        newOne.indices_[newOne.nElements_++] = indexValue;
      } else {
        newOne.elements_[indexValue] = 0.0;
      }
    }
  }
  return newOne;
}
// The sum of two indexed vectors
void CoinIndexedVector::operator+=(const CoinIndexedVector &op2)
{
  // do slowly at first
  *this = *this + op2;
}

// The difference of two indexed vectors
void CoinIndexedVector::operator-=(const CoinIndexedVector &op2)
{
  // do slowly at first
  *this = *this - op2;
}

// The element-wise product of two indexed vectors
void CoinIndexedVector::operator*=(const CoinIndexedVector &op2)
{
  // do slowly at first
  *this = *this * op2;
}

// The element-wise ratio of two indexed vectors (0.0/0.0 => 0.0) (0 vanishes)
void CoinIndexedVector::operator/=(const CoinIndexedVector &op2)
{
  // do slowly at first
  *this = *this / op2;
}
//#############################################################################
void CoinIndexedVector::sortDecrIndex()
{
  // Should replace with std sort
  double *elements = new double[nElements_];
  CoinZeroN(elements, nElements_);
  CoinSort_2(indices_, indices_ + nElements_, elements,
    CoinFirstGreater_2< int, double >());
  delete[] elements;
}

void CoinIndexedVector::sortPacked()
{
  assert(packedMode_);
  CoinSort_2(indices_, indices_ + nElements_, elements_);
}

void CoinIndexedVector::sortIncrElement()
{
  double *elements = new double[nElements_];
  int i;
  for (i = 0; i < nElements_; i++)
    elements[i] = elements_[indices_[i]];
  CoinSort_2(elements, elements + nElements_, indices_,
    CoinFirstLess_2< double, int >());
  delete[] elements;
}

void CoinIndexedVector::sortDecrElement()
{
  double *elements = new double[nElements_];
  int i;
  for (i = 0; i < nElements_; i++)
    elements[i] = elements_[indices_[i]];
  CoinSort_2(elements, elements + nElements_, indices_,
    CoinFirstGreater_2< double, int >());
  delete[] elements;
}
//#############################################################################

void CoinIndexedVector::gutsOfSetVector(int size,
  const int *inds, const double *elems)
{
#ifndef COIN_FAST_CODE
  if (size < 0)
    throw CoinError("negative number of indices", "setVector", "CoinIndexedVector");
#endif
  assert(!packedMode_);
  // find largest
  int i;
  int maxIndex = -1;
  for (i = 0; i < size; i++) {
    int indexValue = inds[i];
#ifndef COIN_FAST_CODE
    if (indexValue < 0)
      throw CoinError("negative index", "setVector", "CoinIndexedVector");
#endif
    if (maxIndex < indexValue)
      maxIndex = indexValue;
  }
  reserve(maxIndex + 1);
  nElements_ = 0;
  // elements_ array is all zero
  bool needClean = false;
  int numberDuplicates = 0;
  for (i = 0; i < size; i++) {
    int indexValue = inds[i];
    if (elements_[indexValue] == 0) {
      if (fabs(elems[i]) >= COIN_INDEXED_TINY_ELEMENT) {
        indices_[nElements_++] = indexValue;
        elements_[indexValue] = elems[i];
      }
    } else {
      numberDuplicates++;
      elements_[indexValue] += elems[i];
      if (fabs(elements_[indexValue]) < COIN_INDEXED_TINY_ELEMENT)
        needClean = true; // need to go through again
    }
  }
  if (needClean) {
    // go through again
    size = nElements_;
    nElements_ = 0;
    for (i = 0; i < size; i++) {
      int indexValue = indices_[i];
      double value = elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        indices_[nElements_++] = indexValue;
      } else {
        elements_[indexValue] = 0.0;
      }
    }
  }
  if (numberDuplicates)
    throw CoinError("duplicate index", "setVector", "CoinIndexedVector");
}

//#############################################################################

void CoinIndexedVector::gutsOfSetVector(int size, int numberIndices,
  const int *inds, const double *elems)
{
  assert(!packedMode_);

  int i;
  reserve(size);
#ifndef COIN_FAST_CODE
  if (numberIndices < 0)
    throw CoinError("negative number of indices", "setVector", "CoinIndexedVector");
#endif
  nElements_ = 0;
  // elements_ array is all zero
  bool needClean = false;
  int numberDuplicates = 0;
  for (i = 0; i < numberIndices; i++) {
    int indexValue = inds[i];
#ifndef COIN_FAST_CODE
    if (indexValue < 0)
      throw CoinError("negative index", "setVector", "CoinIndexedVector");
    else if (indexValue >= size)
      throw CoinError("too large an index", "setVector", "CoinIndexedVector");
#endif
    if (elements_[indexValue]) {
      numberDuplicates++;
      elements_[indexValue] += elems[indexValue];
      if (fabs(elements_[indexValue]) < COIN_INDEXED_TINY_ELEMENT)
        needClean = true; // need to go through again
    } else {
#ifndef COIN_FAC_NEW
      if (fabs(elems[indexValue]) >= COIN_INDEXED_TINY_ELEMENT) {
#endif
        elements_[indexValue] = elems[indexValue];
        indices_[nElements_++] = indexValue;
#ifndef COIN_FAC_NEW
      }
#endif
    }
  }
  if (needClean) {
    // go through again
    size = nElements_;
    nElements_ = 0;
    for (i = 0; i < size; i++) {
      int indexValue = indices_[i];
      double value = elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        indices_[nElements_++] = indexValue;
      } else {
        elements_[indexValue] = 0.0;
      }
    }
  }
  if (numberDuplicates)
    throw CoinError("duplicate index", "setVector", "CoinIndexedVector");
}

//-----------------------------------------------------------------------------

void CoinIndexedVector::gutsOfSetPackedVector(int size, int numberIndices,
  const int *inds, const double *elems)
{
  packedMode_ = true;

  int i;
  reserve(size);
#ifndef COIN_FAST_CODE
  if (numberIndices < 0)
    throw CoinError("negative number of indices", "setVector", "CoinIndexedVector");
#endif
  nElements_ = 0;
  // elements_ array is all zero
  // Does not check for duplicates
  for (i = 0; i < numberIndices; i++) {
    int indexValue = inds[i];
#ifndef COIN_FAST_CODE
    if (indexValue < 0)
      throw CoinError("negative index", "setVector", "CoinIndexedVector");
      //else if (indexValue>=size)
      //throw CoinError("too large an index", "setVector", "CoinIndexedVector");
#endif
    if (fabs(elems[i]) >= COIN_INDEXED_TINY_ELEMENT) {
      elements_[nElements_] = elems[i];
      indices_[nElements_++] = indexValue;
    }
  }
}

//-----------------------------------------------------------------------------

void CoinIndexedVector::gutsOfSetConstant(int size,
  const int *inds, double value)
{

  assert(!packedMode_);
#ifndef COIN_FAST_CODE
  if (size < 0)
    throw CoinError("negative number of indices", "setConstant", "CoinIndexedVector");
#endif
  // find largest
  int i;
  int maxIndex = -1;
  for (i = 0; i < size; i++) {
    int indexValue = inds[i];
#ifndef COIN_FAST_CODE
    if (indexValue < 0)
      throw CoinError("negative index", "setConstant", "CoinIndexedVector");
#endif
    if (maxIndex < indexValue)
      maxIndex = indexValue;
  }

  reserve(maxIndex + 1);
  nElements_ = 0;
  int numberDuplicates = 0;
  // elements_ array is all zero
  bool needClean = false;
  for (i = 0; i < size; i++) {
    int indexValue = inds[i];
    if (elements_[indexValue] == 0) {
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        elements_[indexValue] += value;
        indices_[nElements_++] = indexValue;
      }
    } else {
      numberDuplicates++;
      elements_[indexValue] += value;
      if (fabs(elements_[indexValue]) < COIN_INDEXED_TINY_ELEMENT)
        needClean = true; // need to go through again
    }
  }
  if (needClean) {
    // go through again
    size = nElements_;
    nElements_ = 0;
    for (i = 0; i < size; i++) {
      int indexValue = indices_[i];
      double value = elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        indices_[nElements_++] = indexValue;
      } else {
        elements_[indexValue] = 0.0;
      }
    }
  }
  if (numberDuplicates)
    throw CoinError("duplicate index", "setConstant", "CoinIndexedVector");
}

//#############################################################################
// Append a CoinIndexedVector to the end
void CoinIndexedVector::append(const CoinIndexedVector &caboose)
{
  const int cs = caboose.getNumElements();

  const int *cind = caboose.getIndices();
  const double *celem = caboose.denseVector();
  int maxIndex = -1;
  int i;
  for (i = 0; i < cs; i++) {
    int indexValue = cind[i];
#ifndef COIN_FAST_CODE
    if (indexValue < 0)
      throw CoinError("negative index", "append", "CoinIndexedVector");
#endif
    if (maxIndex < indexValue)
      maxIndex = indexValue;
  }
  reserve(maxIndex + 1);
  bool needClean = false;
  int numberDuplicates = 0;
  for (i = 0; i < cs; i++) {
    int indexValue = cind[i];
    if (elements_[indexValue]) {
      numberDuplicates++;
      elements_[indexValue] += celem[indexValue];
      if (fabs(elements_[indexValue]) < COIN_INDEXED_TINY_ELEMENT)
        needClean = true; // need to go through again
    } else {
      if (fabs(celem[indexValue]) >= COIN_INDEXED_TINY_ELEMENT) {
        elements_[indexValue] = celem[indexValue];
        indices_[nElements_++] = indexValue;
      }
    }
  }
  if (needClean) {
    // go through again
    int size = nElements_;
    nElements_ = 0;
    for (i = 0; i < size; i++) {
      int indexValue = indices_[i];
      double value = elements_[indexValue];
      if (fabs(value) >= COIN_INDEXED_TINY_ELEMENT) {
        indices_[nElements_++] = indexValue;
      } else {
        elements_[indexValue] = 0.0;
      }
    }
  }
  if (numberDuplicates)
    throw CoinError("duplicate index", "append", "CoinIndexedVector");
}
// Append a CoinIndexedVector to the end and modify indices
void CoinIndexedVector::append(CoinIndexedVector &other, int adjustIndex, bool zapElements /*,double multiplier*/)
{
  const int cs = other.nElements_;
  const int *cind = other.indices_;
  double *celem = other.elements_;
  int *newInd = indices_ + nElements_;
  if (packedMode_) {
    double *newEls = elements_ + nElements_;
    if (zapElements) {
      if (other.packedMode_) {
        for (int i = 0; i < cs; i++) {
          newInd[i] = cind[i] + adjustIndex;
          newEls[i] = celem[i] /* *multiplier*/;
          celem[i] = 0.0;
        }
      } else {
        for (int i = 0; i < cs; i++) {
          int k = cind[i];
          newInd[i] = k + adjustIndex;
          newEls[i] = celem[k] /* *multiplier*/;
          celem[k] = 0.0;
        }
      }
    } else {
      if (other.packedMode_) {
        for (int i = 0; i < cs; i++) {
          newEls[i] = celem[i] /* *multiplier*/;
          newInd[i] = cind[i] + adjustIndex;
        }
      } else {
        for (int i = 0; i < cs; i++) {
          int k = cind[i];
          newInd[i] = k + adjustIndex;
          newEls[i] = celem[k] /* *multiplier*/;
        }
      }
    }
  } else {
    double *newEls = elements_ + adjustIndex;
    if (zapElements) {
      if (other.packedMode_) {
        for (int i = 0; i < cs; i++) {
          int k = cind[i];
          newInd[i] = k + adjustIndex;
          newEls[k] = celem[i] /* *multiplier*/;
          celem[i] = 0.0;
        }
      } else {
        for (int i = 0; i < cs; i++) {
          int k = cind[i];
          newInd[i] = k + adjustIndex;
          newEls[k] = celem[k] /* *multiplier*/;
          celem[k] = 0.0;
        }
      }
    } else {
      if (other.packedMode_) {
        for (int i = 0; i < cs; i++) {
          int k = cind[i];
          newInd[i] = k + adjustIndex;
          newEls[k] = celem[i] /* *multiplier*/;
        }
      } else {
        for (int i = 0; i < cs; i++) {
          int k = cind[i];
          newInd[i] = k + adjustIndex;
          newEls[k] = celem[k] /* *multiplier*/;
        }
      }
    }
  }
  nElements_ += cs;
  if (zapElements)
    other.nElements_ = 0;
}
#ifndef CLP_NO_VECTOR
/* Equal. Returns true if vectors have same length and corresponding
   element of each vector is equal. */
bool CoinIndexedVector::operator==(const CoinPackedVectorBase &rhs) const
{
  const int cs = rhs.getNumElements();

  const int *cind = rhs.getIndices();
  const double *celem = rhs.getElements();
  if (nElements_ != cs)
    return false;
  int i;
  bool okay = true;
  for (i = 0; i < cs; i++) {
    int iRow = cind[i];
    if (celem[i] != elements_[iRow]) {
      okay = false;
      break;
    }
  }
  return okay;
}
// Not equal
bool CoinIndexedVector::operator!=(const CoinPackedVectorBase &rhs) const
{
  const int cs = rhs.getNumElements();

  const int *cind = rhs.getIndices();
  const double *celem = rhs.getElements();
  if (nElements_ != cs)
    return true;
  int i;
  bool okay = false;
  for (i = 0; i < cs; i++) {
    int iRow = cind[i];
    if (celem[i] != elements_[iRow]) {
      okay = true;
      break;
    }
  }
  return okay;
}
#endif
// Equal with a tolerance (returns -1 or position of inequality).
int CoinIndexedVector::isApproximatelyEqual(const CoinIndexedVector &rhs, double tolerance) const
{
  CoinIndexedVector tempA(*this);
  CoinIndexedVector tempB(rhs);
  int *cind = tempB.indices_;
  double *celem = tempB.elements_;
  double *elem = tempA.elements_;
  int cs = tempB.nElements_;
  int bad = -1;
  CoinRelFltEq eq(tolerance);
  if (!packedMode_ && !tempB.packedMode_) {
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      if (!eq(celem[iRow], elem[iRow])) {
        bad = iRow;
        break;
      } else {
        celem[iRow] = elem[iRow] = 0.0;
      }
    }
    cs = tempA.nElements_;
    cind = tempA.indices_;
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      if (!eq(celem[iRow], elem[iRow])) {
        bad = iRow;
        break;
      } else {
        celem[iRow] = elem[iRow] = 0.0;
      }
    }
  } else if (packedMode_ && tempB.packedMode_) {
    double *celem2 = rhs.elements_;
    memset(celem, 0, std::min(capacity_, tempB.capacity_) * sizeof(double));
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      celem[iRow] = celem2[i];
    }
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      if (!eq(celem[iRow], elem[i])) {
        bad = iRow;
        break;
      } else {
        celem[iRow] = elem[i] = 0.0;
      }
    }
  } else {
    double *celem2 = elem;
    double *celem3 = celem;
    if (packedMode_) {
      celem2 = celem;
      celem3 = elem;
    }
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      if (!eq(celem2[iRow], celem3[i])) {
        bad = iRow;
        break;
      } else {
        celem2[iRow] = celem3[i] = 0.0;
      }
    }
  }
  if (bad < 0) {
    for (int i = 0; i < tempA.capacity_; i++) {
      if (elem[i]) {
        if (fabs(elem[i]) > tolerance) {
          bad = i;
          break;
        }
      }
    }
    for (int i = 0; i < tempB.capacity_; i++) {
      if (celem[i]) {
        if (fabs(celem[i]) > tolerance) {
          bad = i;
          break;
        }
      }
    }
  }
  return bad;
}
/* Equal. Returns true if vectors have same length and corresponding
   element of each vector is equal. */
bool CoinIndexedVector::operator==(const CoinIndexedVector &rhs) const
{
  const int cs = rhs.nElements_;

  const int *cind = rhs.indices_;
  const double *celem = rhs.elements_;
  if (nElements_ != cs)
    return false;
  bool okay = true;
  CoinRelFltEq eq(1.0e-8);
  if (!packedMode_ && !rhs.packedMode_) {
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      if (!eq(celem[iRow], elements_[iRow])) {
        okay = false;
        break;
      }
    }
  } else if (packedMode_ && rhs.packedMode_) {
    double *temp = new double[std::max(capacity_, rhs.capacity_)];
    memset(temp, 0, std::max(capacity_, rhs.capacity_) * sizeof(double));
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      temp[iRow] = celem[i];
    }
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      if (!eq(temp[iRow], elements_[i])) {
        okay = false;
        break;
      }
    }
    delete[] temp;
  } else {
    const double *celem2 = elements_;
    if (packedMode_) {
      celem2 = celem;
      celem = elements_;
    }
    for (int i = 0; i < cs; i++) {
      int iRow = cind[i];
      if (!eq(celem2[iRow], celem[i])) {
        okay = false;
        break;
      }
    }
  }
  return okay;
}
/// Not equal
bool CoinIndexedVector::operator!=(const CoinIndexedVector &rhs) const
{
  const int cs = rhs.nElements_;

  const int *cind = rhs.indices_;
  const double *celem = rhs.elements_;
  if (nElements_ != cs)
    return true;
  int i;
  bool okay = false;
  for (i = 0; i < cs; i++) {
    int iRow = cind[i];
    if (celem[iRow] != elements_[iRow]) {
      okay = true;
      break;
    }
  }
  return okay;
}
// Get value of maximum index
int CoinIndexedVector::getMaxIndex() const
{
  int maxIndex = -COIN_INT_MAX;
  int i;
  for (i = 0; i < nElements_; i++)
    maxIndex = std::max(maxIndex, indices_[i]);
  return maxIndex;
}
// Get value of minimum index
int CoinIndexedVector::getMinIndex() const
{
  int minIndex = COIN_INT_MAX;
  int i;
  for (i = 0; i < nElements_; i++)
    minIndex = std::min(minIndex, indices_[i]);
  return minIndex;
}
// Scan dense region and set up indices
int CoinIndexedVector::scan()
{
  nElements_ = 0;
  return scan(0, capacity_);
}
// Scan dense region from start to < end and set up indices
int CoinIndexedVector::scan(int start, int end)
{
  assert(!packedMode_);
  end = std::min(end, capacity_);
  start = std::max(start, 0);
  int i;
  int number = 0;
  int *indices = indices_ + nElements_;
  for (i = start; i < end; i++)
    if (elements_[i])
      indices[number++] = i;
  nElements_ += number;
  return number;
}
// Scan dense region and set up indices with tolerance
int CoinIndexedVector::scan(double tolerance)
{
  nElements_ = 0;
#if ABOCA_LITE_FACTORIZATION == 0
  return scan(0, capacity_, tolerance);
#else
  return scan(0, capacity_ & 0x7fffffff, tolerance);
#endif
}
// Scan dense region from start to < end and set up indices with tolerance
int CoinIndexedVector::scan(int start, int end, double tolerance)
{
  assert(!packedMode_);
#if ABOCA_LITE_FACTORIZATION == 0
  end = std::min(end, capacity_);
#else
  end = std::min(end, capacity_ & 0x7fffffff);
#endif
  start = std::max(start, 0);
  int i;
  int number = 0;
  int *indices = indices_ + nElements_;
  for (i = start; i < end; i++) {
    double value = elements_[i];
    if (value) {
      if (fabs(value) >= tolerance)
        indices[number++] = i;
      else
        elements_[i] = 0.0;
    }
  }
  nElements_ += number;
  return number;
}
// These pack down
int CoinIndexedVector::cleanAndPack(double tolerance)
{
  if (!packedMode_) {
    int number = nElements_;
    int i;
    nElements_ = 0;
    for (i = 0; i < number; i++) {
      int indexValue = indices_[i];
      double value = elements_[indexValue];
      elements_[indexValue] = 0.0;
      if (fabs(value) >= tolerance) {
        elements_[nElements_] = value;
        indices_[nElements_++] = indexValue;
      }
    }
    packedMode_ = true;
  }
  return nElements_;
}
// These pack down
int CoinIndexedVector::cleanAndPackSafe(double tolerance)
{
  int number = nElements_;
  if (number) {
    int i;
    nElements_ = 0;
    assert(!packedMode_);
    double *temp = NULL;
    bool gotMemory;
    if (number * 3 < capacity_ - 3 - 9999999) {
      // can find room without new
      gotMemory = false;
      // But may need to align on 8 byte boundary
      char *tempC = reinterpret_cast< char * >(indices_ + number);
      CoinInt64 xx = reinterpret_cast< CoinInt64 >(tempC);
      CoinInt64 iBottom = xx & 7;
      if (iBottom)
        tempC += 8 - iBottom;
      temp = reinterpret_cast< double * >(tempC);
      xx = reinterpret_cast< CoinInt64 >(temp);
      iBottom = xx & 7;
      assert(!iBottom);
    } else {
      // might be better to do complete scan
      gotMemory = true;
      temp = new double[number];
    }
    for (i = 0; i < number; i++) {
      int indexValue = indices_[i];
      double value = elements_[indexValue];
      elements_[indexValue] = 0.0;
      if (fabs(value) >= tolerance) {
        temp[nElements_] = value;
        indices_[nElements_++] = indexValue;
      }
    }
    CoinMemcpyN(temp, nElements_, elements_);
    if (gotMemory)
      delete[] temp;
    packedMode_ = true;
  }
  return nElements_;
}
// Scan dense region and set up indices
int CoinIndexedVector::scanAndPack()
{
  nElements_ = 0;
  return scanAndPack(0, capacity_);
}
// Scan dense region from start to < end and set up indices
int CoinIndexedVector::scanAndPack(int start, int end)
{
  assert(!packedMode_);
  end = std::min(end, capacity_);
  start = std::max(start, 0);
  int i;
  int number = 0;
  int *indices = indices_ + nElements_;
  for (i = start; i < end; i++) {
    double value = elements_[i];
    elements_[i] = 0.0;
    if (value) {
      elements_[number] = value;
      indices[number++] = i;
    }
  }
  nElements_ += number;
  packedMode_ = true;
  return number;
}
// Scan dense region and set up indices with tolerance
int CoinIndexedVector::scanAndPack(double tolerance)
{
  nElements_ = 0;
  return scanAndPack(0, capacity_, tolerance);
}
// Scan dense region from start to < end and set up indices with tolerance
int CoinIndexedVector::scanAndPack(int start, int end, double tolerance)
{
  assert(!packedMode_);
  end = std::min(end, capacity_);
  start = std::max(start, 0);
  int i;
  int number = 0;
  int *indices = indices_ + nElements_;
  for (i = start; i < end; i++) {
    double value = elements_[i];
    elements_[i] = 0.0;
    if (fabs(value) >= tolerance) {
      elements_[number] = value;
      indices[number++] = i;
    }
  }
  nElements_ += number;
  packedMode_ = true;
  return number;
}
// This is mainly for testing - goes from packed to indexed
void CoinIndexedVector::expand()
{
  if (nElements_ && packedMode_) {
    double *temp = new double[capacity_];
    int i;
    for (i = 0; i < nElements_; i++)
      temp[indices_[i]] = elements_[i];
    CoinZeroN(elements_, nElements_);
    for (i = 0; i < nElements_; i++) {
      int iRow = indices_[i];
      elements_[iRow] = temp[iRow];
    }
    delete[] temp;
  }
  packedMode_ = false;
}
// Create packed array
void CoinIndexedVector::createPacked(int number, const int *indices,
  const double *elements)
{
  nElements_ = number;
  packedMode_ = true;
  CoinMemcpyN(indices, number, indices_);
  CoinMemcpyN(elements, number, elements_);
}
// Create packed array
void CoinIndexedVector::createUnpacked(int number, const int *indices,
  const double *elements)
{
  nElements_ = number;
  packedMode_ = false;
  for (int i = 0; i < nElements_; i++) {
    int iRow = indices[i];
    indices_[i] = iRow;
    elements_[iRow] = elements[i];
  }
}
// Create unpacked singleton
void CoinIndexedVector::createOneUnpackedElement(int index, double element)
{
  nElements_ = 1;
  packedMode_ = false;
  indices_[0] = index;
  elements_[index] = element;
}
//  Print out
void CoinIndexedVector::print() const
{
  printf("Vector has %d elements (%spacked mode)\n", nElements_, packedMode_ ? "" : "un");
  for (int i = 0; i < nElements_; i++) {
    if (i && (i % 5 == 0))
      printf("\n");
    int index = indices_[i];
    double value = packedMode_ ? elements_[i] : elements_[index];
    printf(" (%d,%g)", index, value);
  }
  printf("\n");
}

// Zero out array
void CoinArrayWithLength::clear()
{
  assert((size_ > 0 && array_) || !array_);
  memset(array_, 0, size_);
}
// Get array with alignment
void CoinArrayWithLength::getArray(CoinByteArray size)
{
  if (size > 0) {
    if (alignment_ > 2) {
      offset_ = 1 << alignment_;
    } else {
      offset_ = 0;
    }
    assert(size > 0);
    char *array = new char[size + offset_];
    if (offset_) {
      // need offset
      CoinInt64 xx = reinterpret_cast< CoinInt64 >(array);
      int iBottom = static_cast< int >(xx & ((offset_ - 1)));
      if (iBottom)
        offset_ = offset_ - iBottom;
      else
        offset_ = 0;
      array_ = array + offset_;
      ;
    } else {
      array_ = array;
    }
    if (size_ != -1)
      size_ = size;
  } else {
    array_ = NULL;
  }
}
// Get rid of array with alignment
void CoinArrayWithLength::conditionalDelete()
{
  if (size_ == -1) {
    char *charArray = reinterpret_cast< char * >(array_);
    if (charArray)
      delete[](charArray - offset_);
    array_ = NULL;
  } else if (size_ >= 0) {
    size_ = -size_ - 2;
  }
}
// Really get rid of array with alignment
void CoinArrayWithLength::reallyFreeArray()
{
  char *charArray = reinterpret_cast< char * >(array_);
  if (charArray)
    delete[](charArray - offset_);
  array_ = NULL;
  size_ = -1;
}
// Get enough space
void CoinArrayWithLength::getCapacity(CoinByteArray numberBytes, CoinByteArray numberNeeded)
{
  CoinByteArray k = capacity();
  if (k < numberBytes) {
    CoinByteArray saveSize = size_;
    reallyFreeArray();
    size_ = saveSize;
    getArray(std::max(numberBytes, numberNeeded));
  } else if (size_ < 0) {
    size_ = -size_ - 2;
  }
}
/* Alternate Constructor - length in bytes 
   mode -  0 size_ set to size
   >0 size_ set to size and zeroed
   if size<=0 just does alignment
   If abs(mode) >2 then align on that as power of 2
*/
CoinArrayWithLength::CoinArrayWithLength(CoinByteArray size, int mode)
{
  alignment_ = abs(mode);
  size_ = size;
  getArray(size);
  if (mode > 0 && array_)
    memset(array_, 0, size);
}
CoinArrayWithLength::~CoinArrayWithLength()
{
  if (array_)
    delete[](array_ - offset_);
}
// Conditionally gets new array
char *
CoinArrayWithLength::conditionalNew(CoinByteArray sizeWanted)
{
  if (size_ == -1) {
    getCapacity(sizeWanted);
  } else {
    CoinByteArray newSize = (sizeWanted * 101 / 100) + 64;
    // round to multiple of 16
    newSize -= newSize & 15;
    getCapacity(sizeWanted, newSize);
  }
  return array_;
}
/* Copy constructor. */
CoinArrayWithLength::CoinArrayWithLength(const CoinArrayWithLength &rhs)
{
  assert(capacity() >= 0);
  getArray(rhs.capacity());
  if (size_ > 0)
    CoinMemcpyN(rhs.array_, size_, array_);
}

/* Copy constructor.2 */
CoinArrayWithLength::CoinArrayWithLength(const CoinArrayWithLength *rhs)
{
  assert(rhs->capacity() >= 0);
  size_ = rhs->size_;
  getArray(rhs->capacity());
  if (size_ > 0)
    CoinMemcpyN(rhs->array_, size_, array_);
}
/* Assignment operator. */
CoinArrayWithLength &
CoinArrayWithLength::operator=(const CoinArrayWithLength &rhs)
{
  if (this != &rhs) {
    assert(rhs.size_ != -1 || !rhs.array_);
    if (rhs.size_ == -1) {
      reallyFreeArray();
    } else {
      getCapacity(rhs.size_);
      if (size_ > 0)
        CoinMemcpyN(rhs.array_, size_, array_);
    }
  }
  return *this;
}
/* Assignment with length (if -1 use internal length) */
void CoinArrayWithLength::copy(const CoinArrayWithLength &rhs, CoinByteArray numberBytes)
{
  if (numberBytes == -1 || numberBytes <= rhs.capacity()) {
    CoinArrayWithLength::operator=(rhs);
  } else {
    assert(numberBytes >= 0);
    getCapacity(numberBytes);
    if (rhs.array_)
      CoinMemcpyN(rhs.array_, numberBytes, array_);
  }
}
/* Assignment with length - does not copy */
void CoinArrayWithLength::allocate(const CoinArrayWithLength &rhs, CoinByteArray numberBytes)
{
  if (numberBytes == -1 || numberBytes <= rhs.capacity()) {
    assert(rhs.size_ != -1 || !rhs.array_);
    if (rhs.size_ == -1) {
      reallyFreeArray();
    } else {
      getCapacity(rhs.size_);
    }
  } else {
    assert(numberBytes >= 0);
    if (size_ == -1) {
      delete[] array_;
      array_ = NULL;
    } else {
      size_ = -1;
    }
    if (rhs.size_ >= 0)
      size_ = numberBytes;
    assert(numberBytes >= 0);
    assert(!array_);
    if (numberBytes)
      array_ = new char[numberBytes];
  }
}
// Does what is needed to set persistence
void CoinArrayWithLength::setPersistence(int flag, CoinByteArray currentLength)
{
  if (flag) {
    if (size_ == -1) {
      if (currentLength && array_) {
        size_ = currentLength;
      } else {
        conditionalDelete();
        size_ = 0;
        array_ = NULL;
      }
    }
  } else {
    size_ = -1;
  }
}
// Swaps memory between two members
void CoinArrayWithLength::swap(CoinArrayWithLength &other)
{
#ifdef COIN_DEVELOP
  if (!(size_ == other.size_ || size_ == -1 || other.size_ == -1))
    printf("Two arrays have sizes - %d and %d\n", size_, other.size_);
#endif
  assert(alignment_ == other.alignment_);
  char *swapArray = other.array_;
  other.array_ = array_;
  array_ = swapArray;
  CoinByteArray swapSize = other.size_;
  other.size_ = size_;
  size_ = swapSize;
  int swapOffset = other.offset_;
  other.offset_ = offset_;
  offset_ = swapOffset;
}
// Extend a persistent array keeping data (size in bytes)
void CoinArrayWithLength::extend(CoinByteArray newSize)
{
  //assert (newSize>=capacity()&&capacity()>=0);
  assert(size_ >= 0); // not much point otherwise
  if (newSize > size_) {
    char *temp = array_;
    getArray(newSize);
    if (temp) {
      CoinMemcpyN(array_, size_, temp);
      delete[](temp - offset_);
    }
    size_ = newSize;
  }
}

/* Default constructor */
CoinPartitionedVector::CoinPartitionedVector()
  : CoinIndexedVector()
{
  memset(startPartition_, 0, ((&numberPartitions_ - startPartition_) + 1) * sizeof(int));
}
/* Copy constructor. */
CoinPartitionedVector::CoinPartitionedVector(const CoinPartitionedVector &rhs)
  : CoinIndexedVector(rhs)
{
  memcpy(startPartition_, rhs.startPartition_, ((&numberPartitions_ - startPartition_) + 1) * sizeof(int));
}
/* Copy constructor.2 */
CoinPartitionedVector::CoinPartitionedVector(const CoinPartitionedVector *rhs)
  : CoinIndexedVector(rhs)
{
  memcpy(startPartition_, rhs->startPartition_, ((&numberPartitions_ - startPartition_) + 1) * sizeof(int));
}
/* Assignment operator. */
CoinPartitionedVector &CoinPartitionedVector::operator=(const CoinPartitionedVector &rhs)
{
  if (this != &rhs) {
    CoinIndexedVector::operator=(rhs);
    memcpy(startPartition_, rhs.startPartition_, ((&numberPartitions_ - startPartition_) + 1) * sizeof(int));
  }
  return *this;
}
/* Destructor */
CoinPartitionedVector::~CoinPartitionedVector()
{
}
// Add up number of elements in partitions
void CoinPartitionedVector::computeNumberElements()
{
  if (numberPartitions_) {
    assert(packedMode_);
    int n = 0;
    for (int i = 0; i < numberPartitions_; i++)
      n += numberElementsPartition_[i];
    nElements_ = n;
  }
}
// Add up number of elements in partitions and pack and get rid of partitions
void CoinPartitionedVector::compact()
{
  if (numberPartitions_) {
    int n = numberElementsPartition_[0];
    numberElementsPartition_[0] = 0;
    for (int i = 1; i < numberPartitions_; i++) {
      int nThis = numberElementsPartition_[i];
      int start = startPartition_[i];
      memmove(indices_ + n, indices_ + start, nThis * sizeof(int));
      memmove(elements_ + n, elements_ + start, nThis * sizeof(double));
      n += nThis;
    }
    nElements_ = n;
    // clean up
    for (int i = 1; i < numberPartitions_; i++) {
      int nThis = numberElementsPartition_[i];
      int start = startPartition_[i];
      numberElementsPartition_[i] = 0;
      int end = nThis + start;
      if (nElements_ < end) {
        int offset = std::max(nElements_ - start, 0);
        start += offset;
        nThis -= offset;
        memset(elements_ + start, 0, nThis * sizeof(double));
      }
    }
    packedMode_ = true;
    numberPartitions_ = 0;
  }
}
/* Reserve space.
 */
void CoinPartitionedVector::reserve(int n)
{
  CoinIndexedVector::reserve(n);
  memset(startPartition_, 0, ((&numberPartitions_ - startPartition_) + 1) * sizeof(int));
  startPartition_[1] = capacity_; // for safety
}
// Setup partitions (needs end as well)
void CoinPartitionedVector::setPartitions(int number, const int *starts)
{
  if (number) {
    packedMode_ = true;
    assert(number <= COIN_PARTITIONS);
    memcpy(startPartition_, starts, (number + 1) * sizeof(int));
    numberPartitions_ = number;
#ifndef NDEBUG
    assert(startPartition_[0] == 0);
    int last = -1;
    for (int i = 0; i < numberPartitions_; i++) {
      assert(startPartition_[i] >= last);
      assert(numberElementsPartition_[i] == 0);
      last = startPartition_[i];
    }
    assert(startPartition_[numberPartitions_] >= last && startPartition_[numberPartitions_] <= capacity_);
#endif
  } else {
    clearAndReset();
  }
}
// Reset the vector (as if were just created an empty vector). Gets rid of partitions
void CoinPartitionedVector::clearAndReset()
{
  if (numberPartitions_) {
    assert(packedMode_ || !nElements_);
    for (int i = 0; i < numberPartitions_; i++) {
      int n = numberElementsPartition_[i];
      memset(elements_ + startPartition_[i], 0, n * sizeof(double));
      numberElementsPartition_[i] = 0;
    }
  } else {
    memset(elements_, 0, nElements_ * sizeof(double));
  }
  nElements_ = 0;
  numberPartitions_ = 0;
  startPartition_[1] = capacity_;
  packedMode_ = false;
}
// Reset the vector (as if were just created an empty vector). Keeps partitions
void CoinPartitionedVector::clearAndKeep()
{
  assert(packedMode_);
  for (int i = 0; i < numberPartitions_; i++) {
    int n = numberElementsPartition_[i];
    memset(elements_ + startPartition_[i], 0, n * sizeof(double));
    numberElementsPartition_[i] = 0;
  }
  nElements_ = 0;
}
// Clear a partition.
void CoinPartitionedVector::clearPartition(int partition)
{
  assert(packedMode_);
  assert(partition < COIN_PARTITIONS);
  int n = numberElementsPartition_[partition];
  memset(elements_ + startPartition_[partition], 0, n * sizeof(double));
  numberElementsPartition_[partition] = 0;
}
#ifndef NDEBUG
// For debug check vector is clear i.e. no elements
void CoinPartitionedVector::checkClear()
{
#ifndef NDEBUG
  //printf("checkClear %p\n",this);
  assert(!nElements_);
  //assert(!packedMode_);
  int i;
  for (i = 0; i < capacity_; i++) {
    assert(!elements_[i]);
  }
#else
  if (nElements_) {
    printf("%d nElements_ - checkClear\n", nElements_);
    abort();
  }
  int i;
  int n = 0;
  int k = -1;
  for (i = 0; i < capacity_; i++) {
    if (elements_[i]) {
      n++;
      if (k < 0)
        k = i;
    }
  }
  if (n) {
    printf("%d elements, first %d - checkClear\n", n, k);
    abort();
  }
#endif
}
// For debug check vector is clean i.e. elements match indices
void CoinPartitionedVector::checkClean()
{
  //printf("checkClean %p\n",this);
  int i;
  if (!nElements_) {
    checkClear();
  } else {
    if (packedMode_) {
      for (i = 0; i < nElements_; i++)
        assert(elements_[i]);
      for (; i < capacity_; i++)
        assert(!elements_[i]);
    } else {
      double *copy = new double[capacity_];
      CoinMemcpyN(elements_, capacity_, copy);
      for (i = 0; i < nElements_; i++) {
        int indexValue = indices_[i];
        assert(copy[indexValue]);
        copy[indexValue] = 0.0;
      }
      for (i = 0; i < capacity_; i++)
        assert(!copy[i]);
      delete[] copy;
    }
#ifndef NDEBUG
    // check mark array zeroed
    char *mark = reinterpret_cast< char * >(indices_ + capacity_);
    for (i = 0; i < capacity_; i++) {
      assert(!mark[i]);
    }
#endif
  }
}
#endif
// Scan dense region and set up indices (returns number found)
int CoinPartitionedVector::scan(int partition, double tolerance)
{
  assert(packedMode_);
  assert(partition < COIN_PARTITIONS);
  int n = 0;
  int start = startPartition_[partition];
  double *COIN_RESTRICT elements = elements_ + start;
  int *COIN_RESTRICT indices = indices_ + start;
  int sizePartition = startPartition_[partition + 1] - start;
  if (!tolerance) {
    for (int i = 0; i < sizePartition; i++) {
      double value = elements[i];
      if (value) {
        elements[i] = 0.0;
        elements[n] = value;
        indices[n++] = i + start;
      }
    }
  } else {
    for (int i = 0; i < sizePartition; i++) {
      double value = elements[i];
      if (value) {
        elements[i] = 0.0;
        if (fabs(value) > tolerance) {
          elements[n] = value;
          indices[n++] = i + start;
        }
      }
    }
  }
  numberElementsPartition_[partition] = n;
  return n;
}
//  Print out
void CoinPartitionedVector::print() const
{
  printf("Vector has %d elements (%d partitions)\n", nElements_, numberPartitions_);
  if (!numberPartitions_) {
    CoinIndexedVector::print();
    return;
  }
  double *tempElements = CoinCopyOfArray(elements_, capacity_);
  int *tempIndices = CoinCopyOfArray(indices_, capacity_);
  for (int partition = 0; partition < numberPartitions_; partition++) {
    printf("Partition %d has %d elements\n", partition, numberElementsPartition_[partition]);
    int start = startPartition_[partition];
    double *COIN_RESTRICT elements = tempElements + start;
    int *COIN_RESTRICT indices = tempIndices + start;
    CoinSort_2(indices, indices + numberElementsPartition_[partition], elements);
    for (int i = 0; i < numberElementsPartition_[partition]; i++) {
      if (i && (i % 5 == 0))
        printf("\n");
      int index = indices[i];
      double value = elements[i];
      printf(" (%d,%g)", index, value);
    }
    printf("\n");
  }
  delete[] tempElements;
  delete[] tempIndices;
}
/* Sort the indexed storage vector (increasing indices). */
void CoinPartitionedVector::sort()
{
  assert(packedMode_);
  for (int partition = 0; partition < numberPartitions_; partition++) {
    int start = startPartition_[partition];
    double *COIN_RESTRICT elements = elements_ + start;
    int *COIN_RESTRICT indices = indices_ + start;
    CoinSort_2(indices, indices + numberElementsPartition_[partition], elements);
  }
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/
